Multiobjective optimization to assess dengue control costs using a climate-dependent epidemiological model

Arboviruses, diseases transmitted by arthropods, have become a significant challenge for public health managers. The World Health Organization highlights dengue as responsible for millions of infections worldwide annually. As there is no specific treatment for the disease and no free-of-charge vaccine for mass use in Brazil, the best option is the measures to combat the vector, the Aedes aegypti mosquito. Therefore, we proposed an epidemiological model dependent on temperature, precipitation, and humidity, considering symptomatic and asymptomatic dengue infections. Through computer simulations, we aimed to minimize the amount of insecticides and the social cost demanded to treat patients. We proposed a case study in which our model is fitted with real data from symptomatic dengue-infected humans in an epidemic year in a Brazilian city. Our multiobjective optimization model considers an additional control using larvicide, adulticide, and ultra-low volume spraying. The work’s main contribution is studying the monetary cost of the actions to combat the vector demand versus the hospital cost per confirmed infected, comparing approaches with and without additional control. Results showed that the additional vector control measures are cheaper than the hospital treatment without the vector control would be.

www.nature.com/scientificreports/ in which µ, N, ξ , ψ, η, ν, ǫ, θ M , θ I , φ, C, γ , σ , α, β, µ A , µ F , u A , u F ∈ R + and F S + F E + F I = F. The regions of the system Eq. (1) with biological sense are defined by: that is, the populations are nonnegative, the immature stage population A does not exceed the carrying capacity C of the medium, and the human populations do not exceed the total number N of humans. By the next generation matrix method 13 (see Supplementary Information), the expressions found for the basic offspring number, Q 0 , and the basic reproduction number, R 0 , are: Considering the values of the epidemiological parameters of this system invariant in time and assuming the existence of the mosquito population, two equilibrium points can be determined, given by: 1. P 0 = (N, 0, 0, 0, 0,Ā,F S , 0, 0) , which is the trivial or disease-free equilibrium point, Ā = C 1 − 1 Q 0 and F S = σ α µ F +u FĀ . In this case, both the mosquito and human populations are free of the dengue virus. 2. P 1 = (S * * , E * * , M * * , I * * , R * * , A * * , F * * S , F * * E , F * * I ) , which is the nontrivial equilibrium or epidemic equilibrium point, i.e., when the dengue virus is present. In this case, The disease-free equilibrium point, P 0 = (N, 0, 0, 0, 0,Ā,F S , 0, 0) , is locally asymptotically stable for 0 < R 0 < 1 and unstable for R 0 > 1 . On the other hand, the epidemic equilibrium point, P 1 = (S * * , E * * , M * * , I * * , R * * , A * * , F * * S , F * * E , F * * I ) , is locally asymptotically stable for R 0 > 1 and unstable for 0 < R 0 < 1 . Stability analysis demonstrations are also in the Supplementary Information. Furthermore, we use the fourth-order Runge-Kutta method for computational simulation of the dynamical system Eq. (1), with step = 0.01.
Adjusting the climate-dependent parameters. Table 1 contains the interpretation of all model parameters, showing those that vary over time, their acceptable ranges according to the literature, and those constant or climate-independent. Next, we explain how the parameterization was performed according to each climate variable. To visualize the variation of these parameters in time, see Supplementary Fig. S2 online.
Precipitation-dependent entomological parameter. It is known that the influence of rain is difficult to assess, especially during outbreaks of vector-borne diseases 33,34 . On the other hand, many authors agree that precipitation directly influences the aquatic stages of the vector Ae. aegypti, especially in the development of the larvae. They can reach the adult stage because minimal contact with water is required for this development process. However, there is also the phenomenon that the breeding sites can be washed away if there is a high incidence of rain, which eliminates the aquatic phases of the breeding sites and interrupts the reproductive cycle 33,35 . Thus, only the influence of precipitation on mortality in the aquatic phase, µ A , was considered. The parameterization was done linearly, according to the expression below, adapted from 36 .   , I * * = ηνξ ψF * * I µN(ν + µ) (θ I + µ)(ξ ψF * * I + µN) , R * * = θ M M * * + θ I I * * µ , A * * = C − C(α + µ A + u A )(µ F + u F ) σ αǫφ , , and F * * I = γ F * * E N(µ F + u F ) .
in which P = P(t) is the daily precipitation; P 1 (t) and P 2 (t) represent, respectively, the lowest and highest precipitation values during the study period; µ Amin and µ Amax represent, respectively, the minimum and maximum value of µ A , according to data in the literature. In addition to precipitation, µ A is also influenced by temperature, as described below.
Temperature-dependent entomological parameters. In the literature, it is possible to find that temperature influences some of the entomological parameters 29,33,37 . For example, the authors 29 conducted temperature-controlled laboratory experiments to see under what temperature conditions vector populations thrive. Thus, they established a polynomial for each entomological parameter influenced by temperature, minimizing the adjustment error by the linear least squares method. In this work, the influence of temperature is expressed through the function θ(T) , Eq. (5). It represents the polynomial fit made for each temperature-dependent parameter θ j (T) , according to the experiments done by 29 . For the adjustment according to the temperature of each parameter, a polynomial of degree κ is given by: in which θ 1 (T), . . . , θ 4 (T) are the temperature-dependent components of φ , α , µ A , and µ F , respectively; T is the average daily temperature (in • C), and the coefficients b i j , with i = 0, 1, 2, . . . , κ j and j = 1, 2, . . . , 4 , are fitted by the usual linear least squares method. Table 2 shows the values found for each coefficient b i j of the polynomial, with units days −1 ( • C) −i , The degrees κ j are set to 4 for parameters φ, µ A , and µ F and set to 7 for parameter α.
In the specific case of mortality in the aquatic phase, µ A , which is influenced by temperature and precipitation 3,10 , the dependency can assume a well-behaved function, so that µ A (P, T) = a 0 + a 1 µ A (P) + b 1 µ A (T) + O(P 2 , T 2 ) 11 . For simplicity, as a working hypothesis, and because the model (daily) rates are typically less than one, we retain the linear part for the association. It complies with the literature values, respecting the range shown in Table 1. The value is then defined here by the average values obtained in Eqs. (4) and (5).  37 , i.e., the higher the temperature, the higher the rate γ and the shorter the time required to transmit the virus. Thus, Focks et al. 38 showed this relationship with temperature between 12 and 36 • C. Hence, in this work, the following parameterization was used as a function of temperature for γ 31 : Precipitation, temperature, and humidity-dependent entomological parameter. Many interpretations of carrying capacity, C, refer exclusively to the availability of resources in terms of substances necessary for the species to live and reproduce. However, this amount can be understood from a broader point of view as the environmental potential of the species. An example is the enormous capacity of a breeding site full of water and nutrients to produce mosquitoes from larvae at 25 • C. Adult females lay above the water level inside the breeding grounds or where they may be flooded 39,40 . With sufficient rain and favorable humidity conditions, the eggs hatch, generating new larvae 41 .
Thus, precipitation supplies breeding sites and originates new breeding sites in which the larvae develop before becoming adult mosquitoes. However, evaporation tends to reduce these breeding sites. Therefore, due to the strong influence that the amount of water available in breeding grounds and puddles from rain exert on the carrying capacity, also associated with temperature and humidity, it was considered in this work that C varies over time.
How the parameterization of the carrying capacity took place was based on 27 , as defined below: in which C max is the maximum value that all breeding sites can assume. Following the reasoning of 28 , usually, this value corresponds to three times the size of the human population. Like the carrying capacity depends on the amount of water available, H(t), its variation is given as follows: in which H max is the maximum amount of water a breeding site can store without overflowing; �(t) = P(t) − Evap(t) ; P(t) represents the daily precipitation, and Evap(t) is the daily evaporation rate. This rate, in turn, is given by the following expression, adapted from the Ivanov 42 model: in which T is the average daily temperature, Hum is the daily humidity, k represents the Ivanov model constant, and m is a scale factor. This model also considers evapotranspiration, which differs from evaporation because it evaluates the transfer of water to the atmosphere from the soil and plants' transpiration. Although the practical part of measurements related to evapotranspiration is complex 43 , the Ivanov model is commonly used in the literature due to its lower complexity than other methods 27,44,45 . The cosine function represents the epidemic cycle and seasonal influence on dengue cases during the study period.
Multiobjective optimization model. Multiobjective optimization deals with problems that contain at least two objective functions. Usually, these functions conflict with each other in such a way that improving one of them can lead to the degradation of others. As the model Eq. (1) is generic, in the sense that it could be applied to any disease transmitted by the mosquito Ae. aegypti in any city, a trade-off between the cost of vector control could be portrayed by the function J 1 (u A , t A , u F , t F ) versus the hospital cost demanded to treat infected humans, represented in the function J 2 (u A , t A , u F , t F ) . Thus, Eq. (10) represents a general formulation of the multiobjective optimization problem, which is subject to the dynamic system Eq. (1). The problem could also impose other constraints, for example, the minimum and maximum amount of control investment and minimum and maximum time of control application.
Evap(t) = cos(k(25 • C + T(t)) 2 (100 − Hum(t)))m, Table 2. Model temperature-dependent coefficients parameters 29 . www.nature.com/scientificreports/ where s 1 is the relative cost with control in the immature phase of the vector; u A is the intensity of additional larvicides to be applied in the immature phase, starting at time t A until any time t A + τ A ; s 2 the relative cost with control in the adult phase; u F is the intensity of insecticides to be applied in the adult phase, starting at time t F until any time t F + τ F ; s 3 is the hospital cost of infected humans in any time t I ; and I corresponds to the symptomatic infected population, which demands treatment costs. The search for candidate solutions is carried out in the decision variable space in optimization problems. Due to the characteristic of a multiobjective problem being composed of more than one objective function, this process is quite different. When evaluating these functions in the search for candidate solutions, the decision variables will be mapped from the decision variable space to a new search space associated with the objective functions, which are vectors located in the called objective space. This mapping to the objective space is complex and depends on the functions and constraints of the multiobjective problem. Two compromise solutions that are close in the decision variable space will not necessarily be close in the objective space.
Let X be the decision variable space and let F X be the objective space. Let the vector J be described component-wise by J = (J 1 , J 2 ) and let the set F Y represent the image set of region F X . The point x 1 ∈ F X is said to be dominated if there is another point x 2 ∈ F X that is not worse than x 1 for all objectives, as the following relation shows 2 : A feasible solution is said to be nondominated if it is not dominated by any other feasible solution in X. A Pareto-optimal solution is a vector x ∈ F x concerning the objective space. The set of all Pareto-optimal solutions is called Pareto-optimal set. The image of the Pareto-optimal set in the objective space is called Pareto front 2 .
In this work, we use the Non-dominated Sorting Genetic Algorithm II (NSGA-II) 46 as an optimization algorithm to find solutions for the problem. This optimization algorithm classifies and organizes the compromise solutions in nondominated fronts. The solutions in the first front correspond to the best solutions, while the last contains the worst ones. There is no dominance between the solutions located on the same front, but the following fronts contain results that favor less the natural trade-off relationship of the problem. So, it will be necessary to run the optimization algorithm several times to estimate the Pareto-optimal set. The result will be a set of combined nondominated feasible solutions. Table 3 describes the parameters used for the simulations performed with the NSGA-II. The stopping criterion was when the algorithm reached the maximum number of generations.
In summary, there are two main goals in multiobjective optimization: first, to find the set of nondominated solutions in equilibrium with all involved objectives, and second, to find the most diverse solutions possible, ensuring that they are sparsely spaced in the Pareto-optimal region 47 . The solution that will actually be implemented depends on the choice of one or more decision-makers. In this case, public health managers are responsible for decision-making regarding allocating resources for vector control and reducing the risk of exposure of a given population to a virus and the costs arising from probable hospitalization.
In conclusion, the multiobjective approach is equivalent to varying scenarios. Each point on the Pareto front is a trade-off between control interventions and hospital cost, achieved without the need for massive optimization tests to vary control amounts or control application start time. It is up to the decision-maker to choose the point that will best meet its needs concerning the costs involved. In general, for each level of financial resources available, the method finds additional control associated with the respective reduction in the number of infected people and hospital costs.

Study area.
As a case study, we use the model Eq. (1) in the computational simulations and formulate a multiobjective optimization problem to control the vector Ae. aegypti in the city of Belo Horizonte, the capital of the State of Minas Gerais, Brazil. It is located at latitude 19 • 48 ′ 57 ′′ S and longitude 43 • 57 ′ 15 ′′ W, as shown in Fig. 1a 48 . The city has an area of approximately 331 km 2 and is located 716 km away from Brasília, the Federal Capital.
According to the last census data, the population in 2010 was 2,375,151 inhabitants, and the population in 2021 was estimated at 2,530,701 inhabitants 14 . The biome of the city is the Atlantic forest and the humid temperate climate with a dry season (Köppen: Cwb 49 ), that is, with a dry winter, hot, and humid summer. Figure 1b shows the city's series of precipitation, temperature, and humidity during 2019 using the data provided by the Official data on those infected with classic dengue, dengue with signs of alarm, and severe dengue are provided by 51 in epidemiological weeks. Analyzing the historical series of Belo Horizonte since 2010, every 3 years, there has been a dengue epidemic, with the years with the highest number of infected people being 2010 (51,813 cases), 2013 (97,982 cases), 2016 (156,542 cases) and 2019 (116,320 cases). This case study is interested in data on actual dengue infections in Belo Horizonte, especially in epidemic years, to adjust the model with data from the last epidemic that occurred in 2019. However, the vector control actions are daily, and it was decided to follow this time scale for standardization. The distribution of daily cases of the last dengue epidemic in the city is shown in Fig. 2.

Control measures in Belo Horizonte. Following the Brazilian Ministry of Health assumptions, Belo
Horizonte City Hall carried out two types of control applications yearly. First, the focal treatment, which we will call u A , consists of applying larvicide whenever necessary by endemic diseases combat agents. Homes are visited every 2.5 months, as the drug has a residual effect of 2 months, so five cycles of larvicide control application ( u A ) are performed per year. Second, the perifocal treatment, which we will call u F 1 , is an adulticide applied whenever necessary at strategic points, with biweekly visitation and a residual effect of 2 months as well. Moreover, a third way to reduce the number of mosquitoes, although it is not considered a form of vector control, but blocking transmission, is the treatment of ultra-low volume (ULV) spraying that we will call u F 2 (because it also causes a reduction in female population). ULV is applied only when there are many suspected or confirmed dengue cases. The product has virtually no residual effect, as when it is applied, only mosquitoes in flight will be affected.  www.nature.com/scientificreports/ As we will show later, such control actions were insufficient to contain the 2019 epidemic. Therefore, it was decided to simulate in this work an additional control approach to investigate, in the city of Belo Horizonte, possible improvements in the performance of the fight against Ae. aegypti in epidemic years. In practice, the additional control approach consists of applying an additional quantity of kilograms of larvicides ( u A ) and liters of adulticides ( u F = u F 1 + u F 2 ). However, some particularities need to be explained about the additional control approach. First, the larvicide application will be explained, then control via the perifocal application and, shortly after, the ULV spraying.
The proposal is to maintain the periodicity for larvicide control application ( u A ) because the Brazilian Ministry of Health recommends it. So, the duration of each cycle was limited to 73 days, except for the last cycle, which was left with 72 days, totalizing 364 days. In this way, the control cycles follow the same 2.5-month interval established by the Belo Horizonte City Hall. The differential of the proposal is to allow the optimization algorithm to search for the best moment to start applying additional control in each cycle and to define the amount of kilograms demanded.
For perifocal control ( u F 1 ), there is no information on how many cycles are performed per year on average. So, for simplicity, the same number of cycles as the larvicide control was defined for the u F 1 control. That is, both are performed at the same instant of time ( t A = t F 1 ), representing the work of the same endemic diseases combat agent.
It remains to define the proposal for ULV spraying ( u F 2 ). There needs to be more information on how many control cycles are carried out per year, as it is carried out on demand based on an outbreak in the number of dengue cases. So, in this work, the estimated number of mosquitoes will be considered. As expected, the most significant number of mosquitoes is found when infected humans peak. From the model fitting, the application time that coincides with the peak formed between days 50 and 200 was defined for the search space of the optimization algorithm.
The need for additional control. By Eq. (3), R 0 explicitly depends on u F , as it also depends on F S , which in turn is influenced by both additional control parameters, u A and u F . Also consider that u F = u F 1 + u F 2 . Therefore, it can be inferred that R 0 is influenced by vector control actions. Thus, we sought to establish a relationship between R 0 and the model's additional control parameters in such a way as to represent the control intervals that Before presenting this relationship, it is worth noting that climatic variables influence some model parameters over time. However, this variation is not considered to establish the relationship between the controls and the basic reproduction number. Hence, for this step, it was necessary to fix all parameters. So, we computed the climate-dependent parameter's values at each instant of time during the 364 days of the study horizon. Then, we extracted the average of each one of them to bring a more reliable perspective of R 0 because insecticide control is used even when there is no dengue outbreak. Finally, those fixed values were used to establish the relationship shown in Fig. 3.
According to Fig. 3, as there is a red region in which R 0 (u A , u F ) > 1 , it is possible to infer that those actions to combat the vector Ae. aegypti already held in the city were not enough to prevent a dengue epidemic. Therefore, it is necessary to study the combination of controls that will at least guarantee R 0 (u A , u F ) = 1 , in which there is no epidemic. In this sense, Fig. 3 shows that larvicidal control alone is not enough to prevent an epidemic.  www.nature.com/scientificreports/ On the other hand, the application of more adulticide control would be enough to maintain R 0 (u A , u F ) = 1 . This information is essential for defining the amount of the needed additional control actions described below.
Control modeling in Belo Horizonte. We consider two distinct ways to model the control in the numerical simulations, one for the aquatic control u A and the perifocal control u F 1 and another for the ULV spraying u F 2 . Firstly, the descending control was considered for the aquatic phase control and the perifocal control because they have a residual effect of approximately 2 months. This type of descending control is based on the control amount u A in kilograms and the control amount u F 1 in liters. These amounts are reduced at each instant of time, simulating the residual effect of the focal and perifocal treatments' insecticides. So, the control starts at a specific instant of time ( t 0 ) defined by the optimization algorithm, and its effect can still be observed in the following days ( t 0 + τ , with τ = (τ A , τ F 1 ) = 60 ). Therefore, the control action in the last cycle must start no later than day 304 so that all cycles begin and end in the same year. An illustration to exemplify the type of descending control is shown in Fig. 4. In it, u(t) = (u A (t), u F 1 (t)) because, as stated before, the two controls are carried out simultaneously by the same endemic diseases combat agent.
Thus, by the equation of a line, the descending control is formally defined as: In the case of the ULV, the application lasts only 1 day since there is no residual effect on the ULV spraying. Therefore, we considered the step size control represented in Fig. 5. It allows the application of the chemical u F 2 in liters, starting at the instant of time t 0 and ending at the time t 0 + τ (with τ = τ F 2 = 1 ), that is, the duration of the treatment will be limited to 1 day. The optimization algorithm is responsible for choosing the best day for the application. Finally, we can present the multiobjective formulation of the problem of this case study. We consider the trade-off in minimizing the number of females Ae. aegypti with the lowest intensity of control possible so that there are fewer people susceptible to contracting the dengue virus and, at the same time, a lower public hospital cost for those infected. Thus, Eq. (13) represents a general formulation of the multiobjective optimization problem, subject to the constraints of the dynamical system and the bounds of the decision variables in Eq. (14).   www.nature.com/scientificreports/ in which s 1 is the relative cost with control in the immature stage; s 2 the relative cost with perifocal control and s 3 with ULV spraying, both in the adult stage; s 4 the public hospital cost of hospitalizations triggered by classic dengue, dengue with warning signs and severe dengue; u A refers to the intensity of additional larvicides to be applied in the immature phase, starting at time t A ; u F = u F 1 + u F 2 refers to the intensity of additional insecticides to be applied in the adult phase, starting at time t F 1 for perifocal control, and at time t F 2 for ULV spraying; i refers to each control application in each cycle; and I corresponds to the symptomatic infected population, which demands treatment costs. In this work, as the same endemic diseases combat agent applies larvicidal and perifocal control, t A = t F 1 .
The monetary values of costs s 1 , s 2 , s 3 , and s 4 will be described below. All costs mentioned in this work were initially calculated in the local currency (Real -R$) and converted to US dollars (US$) using an average exchange rate (R$ 1 = US$ 0. 25  Before presenting the value of s 4 , which involves the hospital cost, it is necessary to consider the following information. In 2019, 2814 people were hospitalized, 2360 due to classic dengue and 454 due to severe dengue, dengue with warning signs or inconclusive 53 . Among the total infected in the same year, which was 116,320 people, the percentage that required hospitalization in public hospitals was 2.42%. The average cost per hospitalization in public hospitals was US$ 108.78 for classic dengue and US$ 123.22 for dengue with warning signs and severe dengue. Finally, to obtain an average hospital cost per infected person, we calculated: Thus, each infected person cost Belo Horizonte US$ 2.69 in 2019 for public hospital treatment. Across the country, citizens' public hospitalizations are free of charge and funded by the government's Unified Health System (Sistema Único de Saúde-SUS).
Regarding the amount of additional control, following Fig. 3, if the amount of larvicide already applied by the city were doubled, it would still not be enough without the joint application of adulticide (perifocal treatment or ULV spraying). On the other hand, if an additional 20% of adulticide control were applied all the time, for example, even without larvicide, it would guarantee the absence of an epidemic. As it is impracticable to apply control every day, we will allow, following the threshold, up to 20% of an additional control for the control in the immature phase and the one in the adult phase. So, the two controls and the ULV spraying will be applied for an epidemic year. As they are the amount of additional control, in kilograms and liters, u A and u F are intrinsically the sums of the additional control intensity to be used in each cycle of each approach, as the optimization results. (14) subject to:   www.nature.com/scientificreports/ Hence, u A = n 1 u i A kilos and u F = n 1 u i F liters (u F = u F 1 + u F 2 ) , with n representing the number of cycles in each approach and i = 1, . . . , n representing each cycle.
The values the control decision variables can assume in each cycle will now be discussed. In 2019, Belo Horizonte City Hall reported having spent 592 kilos of larvicide, and we restricted up to 20% of this amount for the additional control approach. So, up to 118 kilos were spread over five additional control cycles, totaling approximately 24 kilos or up to 4% per cycle ( 0 u i A 0.04; i = 1, 2, . . . , 5 ). For the adult phase, we divided the 20% additional control into up to 10% perifocal control and up to 10% ULV spraying. In the case of perifocal control, we restricted up to 2% or 0.43 l for additional control in each of the five cycles ( 0 u i F 1 0.02; i = 1, 2, . . . , 5 ). Finally, for ULV spraying, we restricted up to 5 l or 5% for additional control in each of the two cycles ( 0 u i F 2 0.05; i = 1, 2 ). Therefore, the control decision variables' upper bounds in Eq. (14) were chosen to respect the threshold presented in Fig. 3.

Fitting method.
To fit the proposed epidemiological model, we used the lsqcurvefit function of MATLAB® software. The goal is to optimize the recovery rate values of dengue virus for symptomatic ( θ I ) and asymptomatic ( θ M ) infected humans. With this, the evolution of the populations of humans and mosquitoes generated by the numerical simulation of the model Eq. (1) reproduces as faithfully as possible the reality of the populations of humans infected in our case study in the city of Belo Horizonte. Therefore, we also use the set of values of its parameters from Table 1 respecting Belo Horizonte's climate. We used 2019, as the last dengue epidemic occurred, as the base year for control actions.
For the fit to be done, the lsqcurvefit function requires some values and their respective lower and upper limits to be supplied as input parameters. Based on the literature, we set the initial conditions of the parameters to be optimized. In this study, they are the humans' recovery θ I and θ M . The initial conditions of the state variables must also be informed, which were defined here based on information on the number of real cases of infected humans.
Particular attention must be given to the initial condition and limits established for susceptible and recovered humans. From the 116,320 confirmed cases of dengue registered for a population of 2.5 million in 2019, we verified that the incidence of cases was 4.65% of the population. Therefore, the number of 6% of the population, i.e., 150,000 susceptible, was considered reasonable, trying to include also cases of infected people who were not notified, assuming that they did not seek medical help, especially for asymptomatic people. Thus, analyzing the numbers of susceptible and recovered individuals, it can be seen that they make practical sense since the approximately 2.5 million city inhabitants already had the possibility of contact with dengue serotypes for many years. Furthermore, if a serotype infected them in other years, they became immune (recovered) to that serotype.
Information on the relationship between the number of humans and Ae. aegypti mosquitoes is still being determined in the literature. For example, Rodrigues et al. 28 empirically estimated six adult females per human and three larvae per human. In this work, the assumption is that there must be at least one mosquito and up to five for each susceptible. The same was considered for the immature phase of the vector; therefore, the initial condition of 150 thousand was limited to 12.5 million immature and susceptibles adult vector forms. The initial condition was zero for exposed and infected mosquitoes, assuming the absence of dengue virus circulation. In 54 , a survey was conducted where 0.12% of the total number of mosquitoes captured were positive for the dengue virus. So, to also consider the percentage of exposed mosquitoes, we assumed up to 4% of the total susceptibles mosquitoes (i.e., 500,000) as an upper bound for exposed and infected mosquitoes. As for exposed humans, the initial condition is 1000, and the upper bound was considered to be 50,000, dividing the 150,000 of the initial susceptible condition between exposed, symptomatic infected, and asymptomatic infected populations.
Finally, for the fit, we also considered data on precipitation, temperature, and humidity for 2019 in Belo Horizonte. The study horizon contains 364 days, ignoring December 31, based on the 52 epidemiological weeks that a regular year has. Ethical approval. The Belo Horizonte City Hall and the CEFET-MG ethics committees approved this study to access information on how the Ae. aegypti mosquito is controlled in Belo Horizonte. Collecting personal information and participation authorization from city residents was not necessary. Table 4 has the initial condition values and their range, in addition to the optimal values found for the model fit. Figure 6 shows the fit of the simulated symptomatic infected humans to the real curve of confirmed infected. Moreover, Fig. 6 shows the 95% confidence interval for the mean difference between the number of infected estimated by the fitting method and the real infected, which is [− 42.52, 50.95].

Model fitting.
Note how good the fit results were, getting very close to the real curve. Analyzing the fit results for the year 2019, the integral of the curve of the real infected is 116,320, and that of the infected simulated is 121,922, with the root mean squared error equal to 45.4. Figure 7 shows the evolution of human and mosquito populations after curve fitting, representing the final result of the mathematical model. It is essential to mention that in this fit, the control actions carried out in Belo Horizonte are already being considered. However, because of this high number of infected people, this situation gives rise to thinking about the need for additional control, as discussed earlier.
Multiobjective results. After the 30 runs of the NSGA-II, Fig. 8a shows the final Pareto-optimal front. It is composed of 840 combined nondominated points found in all algorithmic runs. Information on the decision variable space and the combined result of all Pareto fronts found in the objective space can be found in the Supplementary Information. www.nature.com/scientificreports/ The Pareto-optimal front has three well-defined regions. First, an almost horizontal region, close to the J 1 function axis; second, an almost vertical region, close to the J 2 function axis; and a third proportional region between the two previous regions, which corresponds to the "knee" of a curve shown in Fig. 8b. Therefore, to enhance comprehension of the obtained combined nondominated front, consider the five highlighted points in Fig. 8.
During the process of choosing these five points, we first select the two extremities points of the Paretooptimal front, that is, where there is the highest hospital cost versus the lowest cost with additional control (point A) and where there is the lowest hospital cost versus the highest cost with additional control (point C). The other points chosen depict options that apparently have a lower hospital cost for an amount of control similar to that of point A and lower cost with control for a hospital cost similar to that of point B. Thus, of the five highlighted points, three (D, B, and E) belong to the knee region highlighted in red, and two (A and C) correspond to the extremes of the other identified regions.
For each of the five chosen points, the control strategies are different. It is up to the decision-maker to choose the strategy best suits its planning. Table 5 shows the percentage of additional control according to each point.
To better understand the costs of control and hospital treatment of those infected at each point, check Table 6, which uses the costs for 2019 as a reference. Consider J 1 + J 2 = J and the efficiency as the reduction of J value.   www.nature.com/scientificreports/ A tiny observation is that the values appearing in the selected points' labels in Fig. 8 are approximations since the software that generates the graphs rounds the values. From the point of view of the efficiency of the objective function value, points E, C and B, in that order, showed the most significant reduction after the additional control compared to the objective function value without additional control. Points B and E are similar in the amount of additional control. One issue that draws attention is the comparison between points A and D. The difference between them is a modest amount of additional perifocal control, showing that this small amount brings a considerable response in the objective function.   www.nature.com/scientificreports/ Given the perifocal control effectiveness, the hospital cost is significantly reduced and impacts the value of J. In this sense, between A and D, the most coherent choice for the decision-maker would be point D, located near the knee region. Figure 9a shows how much could be saved by implementing additional control measures. It also shows reduced hospital costs due to decreased symptomatic infected people. As expected, the total savings are higher for points C and E, at US$ 24,388.72 and US$ 25,018.03, respectively. The difference between the total savings of these points is that in C, more was spent on additional control than in E, but also resulting in more significant hospital cost reduction.
To carry out the sensitivity analysis of the values obtained as total savings, we considered an experiment in which, using the K-nearest neighbors (KNN) method, we found the 10 points closest to each of the five previously selected points (A, B, C, D, and E). Then, we evaluated the decision variables in the dynamic system for each point and obtained the following confidence intervals. Point A-CI: [ Figure 9b shows the percentage of additional control associated with each selected nondominated point. There is a discrepancy regarding the amount of larvicide applied between points C and E. In fact, Fig. 9 shows that the more significant amount found in C further reduces hospital costs but makes the objective function a Nondominated points at the ends of the front (A and C) and the knee point (B). b Nondominated points at the knee region (B, D, and E points).   To conclude the study of the differences between the selected nondominated points, consider Table 7, which reveals the reduction of the symptomatic dengue-infected after additional control actions in each of the five points tested. The efficiency of the reduction of infected was also calculated, consisting of the comparison with the total number of existing infected without additional control.
Concerning the efficiency of the objective function value, points C, E, and B, in that order, showed the most significant reduction after the additional control. In this case, the reduction in the number of infected people who would be hospitalized without additional control actions is more evident at point C, with 156 fewer infected than at point E.
Therefore, the decision maker would probably choose a reduced number of infections among the options presented, which corresponds to the point C control approach. Due to the trade-off nature of the multiobjective problem, it is clear that the reduction of infected is more significant than in point E, although the cost of the additional control operation is slightly higher. However, if there are few financial resources for additional control, the decision maker could choose point E, which would also reduce the number of infected people, but less than point C. Figure 10 shows the estimated temporal dynamics of the mosquito populations after point C interventions to evaluate the impact of this best control measure on the mosquito populations. We added vertical lines to represent the start dates of the implementation of each intervention. So, Fig. 10a shows the five additional larvicide control cycles in the immature phase, the five additional perifocal control cycles, and the two additional ULV spraying cycles in the adult phase of susceptibles mosquitoes. Figure 10b shows the five additional perifocal control cycles and the two additional ULV spraying cycles in the adult phase of exposed and infected mosquitoes. In comparison with Fig. 7d,e, we can see that the most significant reduction was in the exposed and infected mosquitoes population.
The analyzed additional interventions still need to be increased to keep the basic reproduction number below the critical threshold, as shown in Fig. 11 after point C additional control measure implementation. Soon after the end of the residual effect of the control, the system converges to the equilibrium point P 1 , an asymptotically stable node when there is an epidemic to be controlled. The Belo Horizonte City Hall should take more actions to prevent an epidemic, such as a significant number of additional control cycles and effective actions to raise public awareness about removing potential breeding sites. However, the latter is beyond the scope of this work. b Additional control percentage.

Discussion
The nature of the spread of dengue is complex. Understanding it and replicating it in a mathematical model, aiming to obtain predictive information, is a task that involves several fields of knowledge. It is necessary to observe the climatic influence and its variations in precipitation, temperature, and humidity, which as a rule, can change and, consequently, modify the Ae. aegypti life cycle. It is necessary to understand the dynamics between the vector and humans since the mosquito lives in the home environment. It is also necessary to verify whether the decision-making of health managers is being carried out and, if so, with what quality, or even find out if new control measures can or should, be carried out. Many other variables could be mentioned here to reinforce the idea that solving real-world problems is challenging. However, through multiobjective optimization, this work promoted the simultaneous evaluation of some complex variables. Thus, we were able to promote a useful mathematical and computational framework for health managers to the point of generating simplified results and information. For this, we proposed a new epidemiological mathematical model-with emphasis on incorporating asymptomatic populations, which accelerate the dissemination of the dengue virus, and on the aquatic phase of the vector-so that it would be possible to simulate the application of larvicides, as it occurs by the recommendations of the Ministry of Health to face the vector. Some model parameters are influenced by climatic variables directly related to the mosquito life cycle. The model, being generic, can be applied to any city. In this work, we carried State variables 10 5 Fe(t) Fi(t) b Exposed and infected Aedes aegypti mosquitoes populations. www.nature.com/scientificreports/ out a case study for the city of Belo Horizonte, Minas Gerais, Brazil, and therefore, we used the model with real data to represent the curve of infected people in the city during the epidemic years. Next, we presented the multiobjective optimization results using the NSGA-II genetic algorithm. A significant contribution of this work to the literature and society is to have presented the trade-off between the costs involving additional vector control versus the reduced need for public hospital treatment. The additional control approach we proposed was cost-effective from a financial point of view; however, the most relevant finding was the decrease in infected people and the need for hospital treatment compared to the approach without additional control.
In this sense, the health manager will be equipped with information for better decision-making. From now on, lessons learned from the epidemic years can be used to prepare control actions and support awareness campaigns about the risks of dengue (both in the media and in schools). Thus, the population will gain in terms of health, especially the less favored. The lower the need for hospitalization, the lower the exposure to biological risks that can lead to complications in the clinical condition of patients.
During other epidemic years, the health manager can monitor the number of identified dengue cases and then use our methodology to adjust the model to this new scenario. Depending on their needs, the health manager can use smaller or larger time windows than the 1-year period we tested for the 2019 case study. In this way, adopting the estimated optimal control measures and knowing the best moment for implementing the interventions is possible. There is also the mosquito abundance indicator Larval Index Rapid Assay for Aedes aegypti (LIRAa) for entomological surveillance of Aedes, widely used in Brazil by the Ministry of Health. Monitoring this indicator could be coupled in future works to act more effectively in control actions to avoid the beginning of an epidemic.
Limitations of the case study include underreporting of symptomatic infected data and the need for more information on the ideal number of additional control cycles, especially for ULV spraying. When evaluating the effort to apply additional control by Belo Horizonte City Hall, some points were not taken into account, e.g., personnel costs (hiring or overtime and uniforms for endemic diseases combat agents), as well as the cost of machinery, all covered by the City Hall.
There will also be productivity costs associated with people missing paid and unpaid work when they are ill, in treatment, or in premature death. We did not evaluate it in this work, but the additional vector control might save the government more money than just reduced public hospital admission. In Brazil, citizens' public hospitalizations are free of charge and funded by the government. So, another possibility would be to assess dengue's total economic burden by also including information on private hospitalizations, which in Brazil can be funded by companies when they pay for a private health insurance for their employees or by citizens themselves.
Hence, there is no doubt that this work opens up possibilities for solving future problems. This case study only considers insecticide vector control; however, in the future the model could be adapted to include other interventions against the spread of dengue fever, such as Wolbachia bacteria or vaccination. The effect of delayed oviposition and sterile males could also be included. Within the scope of the control application, other approaches could be tested and compared. As dengue is a multifactorial disease, control approaches must be multiple and ecologically structured. In this sense, it is worth exploring the possibility of including the risk of systematic insecticide use and possible environmental impacts in the model. Complementary measures to combat the vector, e.g., intensifying campaigns to raise awareness of the population on mechanical control would also be appropriate. Finally, another interesting approach would be to propose computational models based on cellular automata.
In summary, we develop in this work a novel epidemiological mathematical model and a solution to a multiobjective problem for which we show a set of compromise solutions between the costs of controlling the Ae. aegypti and the number of dengue-infected people considering real data in the case study for the city of Belo Horizonte, which decision-makers can use as a range of options to implement a successful control action according to the financial and social condition of the case study city.